﻿using System;
using System.IO; 

namespace ASquared.ModelOptimization
{
    // Method enum
    public enum RootFinderMethod { Step, Bisection, FalsePosition, Newtons, Secant, Mullers, InvInterpolation, Brents }

    // Main classes
    #region RootFinder interface

    /// <summary>An implementation of a class that will find the roots of a function.</summary>
    public interface IRootFinder
    {
        // Properties
        /// <summary>Gets and sets the function for which to find the root.</summary>
        ISingleVariableFunction Function { get; set; }
        /// <summary>Gets and sets the current inputs for the function for which to calculate the root.</summary>
        RootFinderInputs Inputs { get; set; }
        /// <summary>Gets and sets the current estimate for the root of the function.</summary>
        Double RootEstimate { get; set; }
        /// <summary>Gets and sets he current function value at the current root estimate.</summary>
        Double FunctionValue { get; set; }
        /// <summary>Gets and sets the bounds for the root estimates. By default, these bounds are strictly greater than or less then. If you want "or equal to" change the Strict...Bound properties.</summary>
        Bounds Bounds { get; set;  }
        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        bool IsSolvable { get; }
        /// <summary>Gets a value that specifies the number of iterations that have occurred to estimate the root of the current function.</summary>
        Int32 NumOfIters { get; }

        // Methods [the order goes: set RootEstimate with EstimateRoot(), calculate FunctionValue, check IsConverged(), exit or UpdateInfo()]
        /// <summary>Gets the root of <c>function</c>.</summary>
        /// <returns>Returns the root of <c>function</c>.</returns>
        Double GetRoot();
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        Double EstimateRoot();
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        /// <remarks>EstimateRoot() must be called before this method.</remarks>
        bool IsConverged();
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        /// <remarks>EstimateRoot() must be called before this method.</remarks>
        void UpdateInfo();
        /// <summary>Performs a deep copy of this instance except for the ISingleVariableFunction.</summary>
        IRootFinder Copy();
    }

    #endregion
    #region RootFinderInputs

    /// <summary>Aids in building the inputs array for specific methods.</summary>
    public class RootFinderInputs
    {
        public Double step = Double.NaN, derivVal = Double.NaN,
                    rootInit = Double.NaN, fxnValInit = Double.NaN,
                    rootPrev1 = Double.NaN, fxnValPrev1 = Double.NaN,
                    rootPrev2 = Double.NaN, fxnValPrev2 = Double.NaN;
        
        /// <summary>Copies this instance.</summary>
        public RootFinderInputs Copy()
        {
            return (RootFinderInputs)this.MemberwiseClone();
        }
    }

    #endregion
    #region RootFinderBase

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using a simple stepping method.</summary>
    public class RootFinderBase : IRootFinder
    {
        #region Constructor and log file

        private string logfile = string.Empty;
        private StreamWriter logStream = null; 
        public void WriteLog(string message) 
        {
            if (this.logStream != null)
            {
                this.logStream.WriteLine(message);
                this.logStream.Flush(); 
            }
        }
        public RootFinderBase()
        {
            if (string.IsNullOrEmpty(logfile)) 
            {
                this.logfile = Path.GetDirectoryName(Path.GetFullPath(System.Reflection.Assembly.GetExecutingAssembly().Location)) + @"\rf_log.csv";
                if (File.Exists(this.logfile))
                {
                    try
                    {
                        Stream str = new FileStream(this.logfile, FileMode.OpenOrCreate, FileAccess.Write, FileShare.Write); 
                        this.logStream = new StreamWriter(str);
                    }
                    catch 
                    {
                        this.logStream = null; 
                    }
                }
            }
        }
        ~RootFinderBase()
        {
            if (this.logStream != null)
            {
                try { this.logStream.Close(); }
                catch { }
                this.logStream = null;
            }
        }

        #endregion
        #region Convergence variables and properties

        private Double tolerance = 0.001;
        private Int32 maxIterations = 50;
        /// <summary>Gets and sets the tolerance value for numerical solutions.</summary>
        public Double Tolerance
        {
            get
            {
                return tolerance;
            }
            set
            {
                tolerance = value;
            }
        }
        /// <summary>Gets and sets the maximum number of iterations allowed.</summary>
        public Int32 MaxIterations
        {
            get
            {
                return maxIterations;
            }
            set
            {
                maxIterations = value;
            }
        }

        #endregion
        #region Interface variables and properties

        // Instance variables
        private ISingleVariableFunction fxn = null;
        private Bounds privBounds = new Bounds();
        private RootFinderInputs inputs;
        private Double root, fxnVal;
        private Int32 noOfIters = 0;

        // Properties
        /// <summary>Gets and sets the function for which to calculate the root.</summary>
        public virtual ISingleVariableFunction Function
        {
            get
            {
                return this.fxn;
            }
            set
            {
                this.fxn = value;
            }
        }
        /// <summary>Gets and sets the current estimate for the root of the function.</summary>
        public virtual Double RootEstimate
        {
            get { return this.root; }
            set
            {
                this.root = value;
                this.noOfIters++;
                if (this.noOfIters > maxIterations)
                    throw new Exception("The root-finding algorithm stopped because the maximum number of iterations was reached.");
                this.root = this.Bounds.ClosestBoundedValue(this.root);
                this.WriteLog("root," + this.noOfIters.ToString() + "," + this.root.ToString());
            }
        }
        /// <summary>Gets and sets the current function value evaluated at the current root estimate.</summary>
        public virtual Double FunctionValue
        {
            get { return this.fxnVal; }
            set
            {
                this.fxnVal = value;
                this.WriteLog("fxnVal," + this.noOfIters.ToString() + "," + this.fxnVal); 
            }
        }
        /// <summary>Gets and sets the inputs required for the root-finding algorithm.</summary>
        public virtual RootFinderInputs Inputs
        {
            get
            {
                return this.inputs;
            }
            set
            {
                this.inputs = value;
                this.WriteLog("rootInit," + this.noOfIters.ToString() + "," + this.inputs.rootInit.ToString());
                this.WriteLog("rootPrev1," + this.noOfIters.ToString() + "," + this.inputs.rootPrev1.ToString());
                this.WriteLog("rootPrev2," + this.noOfIters.ToString() + "," + this.inputs.rootPrev2.ToString());
                this.WriteLog("fxnValInit," + this.noOfIters.ToString() + "," + this.inputs.fxnValInit.ToString());
                this.WriteLog("fxnValPrev1," + this.noOfIters.ToString() + "," + this.inputs.fxnValPrev1.ToString());
                this.WriteLog("fxnValPrev2," + this.noOfIters.ToString() + "," + this.inputs.fxnValPrev2.ToString());
                this.WriteLog("derivVal," + this.noOfIters.ToString() + "," + this.inputs.derivVal.ToString());
                this.WriteLog("step," + this.noOfIters.ToString() + "," + this.inputs.step.ToString());
            }
        }
        /// <summary>Gets and sets the bounds on the root estimates. By default, these bounds are strictly less than and greater than. If you want "or equal to", change the <c>StrictBounds</c> property.</summary>
        public virtual Bounds Bounds
        {
            get { return this.privBounds; }
            set { this.privBounds = value; }
        }
        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        public virtual bool IsSolvable
        {
            get { return true; }
        }
        /// <summary>Gets a value that specifies the number of iterations that have occurred to estimate the root of the current function.</summary>
        public virtual Int32 NumOfIters
        {
            get { return this.noOfIters; }
        }

        #endregion
        #region Interface methods

        // Local methods
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public virtual Double GetRoot()
        {
            for (Int32 i = 0; i < maxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = fxn.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public virtual Double EstimateRoot()
        {
            return this.RootEstimate + this.Inputs.step;
        }
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        public virtual bool IsConverged()
        {
            return (Math.Abs(this.FunctionValue) < Tolerance);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public virtual void UpdateInfo()
        {
            if (Math.Abs(this.FunctionValue) > Math.Abs(this.Inputs.fxnValPrev1))
                this.Inputs.step /= -2.0;
            if ((this.Inputs.step < 0 && this.Bounds.IsOnLowerBound(this.RootEstimate)) || (this.Inputs.step > 0 && this.Bounds.IsOnUpperBound(this.RootEstimate)))
                this.Inputs.step = -this.Inputs.step;
            this.Inputs.fxnValPrev1 = this.FunctionValue;
        }
        /// <summary>Performs a deep copy of this instance.</summary>
        public virtual IRootFinder Copy()
        {
            RootFinderBase rf = (RootFinderBase)this.MemberwiseClone();
            rf.Bounds = this.Bounds.Copy();
            rf.Inputs = this.Inputs.Copy();
            return rf;
        }

        #endregion
        #region Static methods

        /// <summary>Gets a root-finder method from a string.</summary>
        /// <param name="rfmethod">The root-finder method as a string.</param>
        public static RootFinderMethod GetMethod(string rfmethod)
        {
            try { return (RootFinderMethod)Enum.Parse(typeof(RootFinderMethod), rfmethod, true); }
            catch { return RootFinderMethod.Secant; } 
        }
        /// <summary>Gets an array of all root-finding methods.</summary>
        public static RootFinderMethod[] GetMethods()
        {
            return (RootFinderMethod[])Enum.GetValues(typeof(RootFinderMethod)); 
        }
        /// <summary>Gets an array of type names for all root-finding methods.</summary>
        public static string[] GetMethodNames()
        {
            return Array.ConvertAll(GetMethods(), method => method.ToString()); 
        }

        #endregion
    }

    #endregion
    #region RootFinder

    /// <summary>An overarching class that will find a root given a particular method.</summary>
    public class RootFinder : RootFinderBase, IRootFinder
    {
        // Static
        #region Static methods

        /// <summary>Retrieves the root-finder class associated with a particular method. If the method cannot be solved with the specified inputs, a different method will automatically be selected (and the new method is given in the output variable: <c>newmethod</c>).</summary>
        /// <param name="method">The desired root-finding method.</param>
        public static IRootFinder GetRootFinder(RootFinderMethod method)
        {
            switch (method)
            {
                case RootFinderMethod.Step: return new RootFinderBase();
                case RootFinderMethod.Bisection: return new BisectionRootFinder();
                case RootFinderMethod.FalsePosition: return new FalsePositionRootFinder();
                case RootFinderMethod.Newtons: return new NewtonsRootFinder();
                case RootFinderMethod.Secant: return new SecantRootFinder();
                case RootFinderMethod.Mullers: return new MullersRootFinder();
                case RootFinderMethod.InvInterpolation: return new InvInterpolationRootFinder();
                case RootFinderMethod.Brents: return new BrentsRootFinder();
                default: throw new Exception("The specified root finder method (\"" + method.ToString() + "\") is undefined.");
            }
        }
        /// <summary>Retrieves the root-finder class associated with a particular method. If the method cannot be solved with the specified inputs, a different method will automatically be selected (and the new method is given in the output variable: <c>newmethod</c>).</summary>
        /// <param name="method">The desired root-finding method.</param>
        /// <param name="Inputs">The inputs for the root-finder class to check if it is solvable.</param>
        /// <param name="newmethod">The method of the returned root-finder class, which may be different than <c>method</c> if it can't be solved given <c>Inputs</c>.</param>
        public static IRootFinder GetRootFinder(RootFinderMethod method, RootFinderInputs Inputs, out RootFinderMethod newmethod)
        {
            IRootFinder rf = GetRootFinder(method);
            rf.Inputs = Inputs;
            newmethod = method;
            if (!rf.IsSolvable)
            {
                newmethod = RootFinderMethod.Secant;
                return new SecantRootFinder(); // use the secant method.
            }
            else
                return rf;
        }

        #endregion

        // Local
        #region Local variables and properties

        private RootFinderMethod method = RootFinderMethod.Brents;
        private IRootFinder rootFinder = new BisectionRootFinder();
        /// <summary>Gets and sets the method to use when updating root value estimates and checking convergence.</summary>
        public RootFinderMethod Method
        {
            get
            {
                return this.method;
            }
            set
            {
                this.method = value;
                this.rootFinder = GetRootFinder(value);
            }
        }

        #endregion

        // Constructor 
        #region Constructors

        /// <summary>Constructs a root-finder based on the specified method.</summary>
        /// <param name="method">Specifies the root-finding method to use.</param>
        /// <param name="InitialInputs">The initial inputs needed for the root-finding method to function.</param>
        public RootFinder(RootFinderMethod method, RootFinderInputs InitialInputs)
            : this(method, InitialInputs, null)
        {
        }
        /// <summary>Constructs a root-finder based on the specified method.</summary>
        /// <param name="method">Specifies the root-finding method to use.</param>
        /// <param name="InitialInputs">The initial inputs needed for the root-finding method to function.</param>
        /// <param name="fxn">The function for which to get the root value.</param>
        public RootFinder(RootFinderMethod method, RootFinderInputs InitialInputs, ISingleVariableFunction fxn)
        {
            this.method = method;
            this.rootFinder = GetRootFinder(method, InitialInputs, out this.method);
            this.Inputs = InitialInputs;
            this.Function = fxn;
        }

        #endregion

        // Interface 
        #region Interface variables and properties

        /// <summary>Gets and sets the function for which to find the root.</summary>
        public override ISingleVariableFunction Function
        {
            get
            {
                return this.rootFinder.Function;
            }
            set
            {
                this.rootFinder.Function = value;
            }
        }
        /// <summary>Gets and sets the estimate of the root value of the function.</summary>
        public override Double RootEstimate
        {
            get
            {
                return this.rootFinder.RootEstimate;
            }
            set
            {
                this.rootFinder.RootEstimate = value;
            }
        }
        /// <summary>Gets and sets the function value at the current root value estimate.</summary>
        public override Double FunctionValue
        {
            get
            {
                return this.rootFinder.FunctionValue;
            }
            set
            {
                this.rootFinder.FunctionValue = value;
            }
        }
        /// <summary>Gets and sets the inputs needed to evaluate the function.</summary>
        public override RootFinderInputs Inputs
        {
            get
            {
                return this.rootFinder.Inputs;
            }
            set
            {
                this.rootFinder.Inputs = value;
            }
        }
        /// <summary>Gets and sets the bounds for the root estimates. By default, these bounds are strictly greater than or less then. If you want "or equal to" change the Strict...Bound properties.</summary>
        public override Bounds Bounds
        {
            get { return this.rootFinder.Bounds; }
            set { this.rootFinder.Bounds = value; }
        }
        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        public override bool IsSolvable
        {
            get
            {
                return this.rootFinder.IsSolvable;
            }
        }
        /// <summary>Gets a value that specifies the number of iterations that have occurred to estimate the root of the current function.</summary>
        public override Int32 NumOfIters
        {
            get
            {
                return this.rootFinder.NumOfIters;
            }
        }

        #endregion
        #region Interface methods

        // Local methods
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            return this.rootFinder.GetRoot();
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            return this.rootFinder.EstimateRoot();
        }
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        public override bool IsConverged()
        {
            return this.rootFinder.IsConverged();
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            this.rootFinder.UpdateInfo();
        }
        /// <summary>Performs a deep copy of this instance except for the ISingleVariableFunction.</summary>
        public override IRootFinder Copy()
        {
            RootFinder rf = (RootFinder)this.MemberwiseClone();
            rf.rootFinder = this.rootFinder.Copy();
            return rf; 
        }

        #endregion
    }

    #endregion

    // Methods
    #region Bisection Method

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the Bisection Method.</summary>
    public class BisectionRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface variables and properties

        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        public override bool IsSolvable
        {
            get
            {
                return !(this.Inputs.fxnValInit * this.Inputs.fxnValPrev1 > 0);
            }
        }

        #endregion
        #region Interface methods

        // Local methods
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            if (!this.IsSolvable)
                throw new Exception("This root-finding method requires an interval that contains the function root value.");
            if (this.Inputs.fxnValInit == 0.0) return this.Inputs.rootInit;
            if (this.Inputs.fxnValPrev1 == 0.0) return this.Inputs.rootPrev1;

            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            return (this.Inputs.rootInit + this.Inputs.rootPrev1) / 2.0;
        }
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        public override bool IsConverged()
        {
            return (Math.Abs(this.FunctionValue) < this.Tolerance || Math.Abs((this.Inputs.rootPrev1 - this.Inputs.rootInit) / 2.0) < this.Tolerance);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            if (this.FunctionValue * this.Inputs.fxnValInit >= 0)
            {
                this.Inputs.rootInit = this.RootEstimate;
                this.Inputs.fxnValInit = this.FunctionValue;
            }
            else
            {
                this.Inputs.rootPrev1 = this.RootEstimate;
                this.Inputs.fxnValPrev1 = this.FunctionValue;
            }
        }

        #endregion
    }

    #endregion
    #region FalsePosition Method

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the False Position Method.</summary>
    public class FalsePositionRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface variables and properties

        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        public override bool IsSolvable
        {
            get
            {
                return !(this.Inputs.fxnValInit * this.Inputs.fxnValPrev1 > 0);
            }
        }

        #endregion
        #region Interface methods

        // Local methods
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            if (!this.IsSolvable)
                throw new Exception("This root-finding method requires an interval that contains the function root value.");
            if (this.Inputs.fxnValInit == 0.0) return this.Inputs.rootInit;
            if (this.Inputs.fxnValPrev1 == 0.0) return this.Inputs.rootPrev1;

            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            return (this.Inputs.rootInit * this.Inputs.fxnValPrev1 - this.Inputs.rootPrev1 * this.Inputs.fxnValInit) / (this.Inputs.fxnValPrev1 - this.Inputs.fxnValInit); 
        }
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        public override bool IsConverged()
        {
            return (Math.Abs(this.FunctionValue) < this.Tolerance || Math.Abs((this.Inputs.rootPrev1 - this.Inputs.rootInit) / 2.0) < this.Tolerance);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            if (this.FunctionValue * this.Inputs.fxnValInit >= 0)
            {
                this.Inputs.rootInit = this.RootEstimate;
                this.Inputs.fxnValInit = this.FunctionValue;
            }
            else
            {
                this.Inputs.rootPrev1 = this.RootEstimate;
                this.Inputs.fxnValPrev1 = this.FunctionValue;
            }
        }

        #endregion
    }

    #endregion
    #region NewtonsMethod

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using Newton's Method.</summary>
    public class NewtonsRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface methods

        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            return this.Inputs.rootInit - this.Inputs.fxnValInit / (this.Inputs.derivVal == 0.0 ? 0.1 : this.Inputs.derivVal);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            ISingleVariableFunction deriv;
            if (this.Function != null && (deriv = this.Function.GetDeriv()) != null)
                this.Inputs.derivVal = deriv.GetValue(this.RootEstimate);
            else
                this.Inputs.derivVal = (this.Inputs.fxnValInit - this.FunctionValue) / (this.Inputs.rootInit - this.RootEstimate);
            this.Inputs.rootInit = this.RootEstimate;
            this.Inputs.fxnValInit = this.FunctionValue;
        }

        #endregion
    }

    #endregion
    #region SecantMethod

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the Secant Method.</summary>
    public class SecantRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface methods

        // Local methods
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            if (this.Inputs.fxnValInit == 0) return this.Inputs.rootInit;
            if (this.Inputs.fxnValPrev1 == 0) return this.Inputs.rootPrev1;
            if (this.Inputs.fxnValInit == this.Inputs.fxnValPrev1) return this.Inputs.rootInit - this.Inputs.fxnValInit;
            return this.Inputs.rootInit - this.Inputs.fxnValInit * (this.Inputs.rootInit - this.Inputs.rootPrev1) / (this.Inputs.fxnValInit - this.Inputs.fxnValPrev1);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            this.Inputs.rootPrev1 = this.Inputs.rootInit;
            this.Inputs.fxnValPrev1 = this.Inputs.fxnValInit;
            this.Inputs.rootInit = this.RootEstimate;
            this.Inputs.fxnValInit = this.FunctionValue;
        }

        #endregion
    }

    #endregion
    #region MullersMethod

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the Secant Method.</summary>
    public class MullersRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface methods

        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            Double w = DividedDifference.Calculate(new Double[] { this.Inputs.rootInit, this.Inputs.rootPrev1 }, new Double[] { this.Inputs.fxnValInit, this.Inputs.fxnValPrev1 })
                + DividedDifference.Calculate(new Double[] { this.Inputs.rootInit, this.Inputs.rootPrev2 }, new Double[] { this.Inputs.fxnValInit, this.Inputs.fxnValPrev2 })
                - DividedDifference.Calculate(new Double[] { this.Inputs.rootPrev1, this.Inputs.rootPrev2 }, new Double[] { this.Inputs.fxnValPrev1, this.Inputs.fxnValPrev2 });
            Double fdiv = DividedDifference.Calculate(new Double[] { this.Inputs.rootInit, this.Inputs.rootPrev1, this.Inputs.rootPrev2 }, new Double[] { this.Inputs.fxnValInit, this.Inputs.fxnValPrev1, this.Inputs.fxnValPrev2 });
            Double inside = Math.Pow(w, 2) - 4 * this.Inputs.fxnValInit * fdiv;
            if (inside < 0) // this will result in unsolvable state... use secant method instead...
                return this.Inputs.rootInit - this.Inputs.fxnValInit * (this.Inputs.rootInit - this.Inputs.rootPrev1) / (this.Inputs.fxnValInit - this.Inputs.fxnValPrev1);
            Double denom = Math.Sqrt(inside);
            if (Math.Abs(w + denom) < Math.Abs(w - denom))
                denom = w - denom;
            else
                denom = w + denom;
            return this.Inputs.rootInit - (2 * this.Inputs.fxnValInit) / denom;
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            this.Inputs.rootPrev2 = this.Inputs.rootPrev1;
            this.Inputs.fxnValPrev2 = this.Inputs.fxnValPrev1;
            this.Inputs.rootPrev1 = this.Inputs.rootInit;
            this.Inputs.fxnValPrev1 = this.Inputs.fxnValInit;
            this.Inputs.rootInit = this.RootEstimate;
            this.Inputs.fxnValInit = this.FunctionValue;
        }

        #endregion
    }

    #endregion
    #region InvInterpolationMethod

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the Secant Method.</summary>
    public class InvInterpolationRootFinder : RootFinderBase, IRootFinder
    {
        #region Interface methods

        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            return this.Inputs.rootPrev2 * this.Inputs.fxnValPrev1 * this.Inputs.fxnValInit / ((this.Inputs.fxnValPrev2 - this.Inputs.fxnValPrev1) * (this.Inputs.fxnValPrev2 - this.Inputs.fxnValInit))
                + this.Inputs.rootPrev1 * this.Inputs.fxnValPrev2 * this.Inputs.fxnValInit / ((this.Inputs.fxnValPrev1 - this.Inputs.fxnValPrev2) * (this.Inputs.fxnValPrev1 - this.Inputs.fxnValInit))
                + this.Inputs.rootInit * this.Inputs.fxnValPrev1 * this.Inputs.fxnValPrev2 / ((this.Inputs.fxnValInit - this.Inputs.fxnValPrev1) * (this.Inputs.fxnValInit - this.Inputs.fxnValPrev2));
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            this.Inputs.rootPrev2 = this.Inputs.rootPrev1;
            this.Inputs.fxnValPrev2 = this.Inputs.fxnValPrev1;
            this.Inputs.rootPrev1 = this.Inputs.rootInit;
            this.Inputs.fxnValPrev1 = this.Inputs.fxnValInit;
            this.Inputs.rootInit = this.RootEstimate;
            this.Inputs.fxnValInit = this.FunctionValue;
        }

        #endregion
    }

    #endregion
    #region BrentsMethod

    /// <summary>A class implementing the <c>IRootFinder</c> interface to find roots of equations using the Secant Method.</summary>
    public class BrentsRootFinder : RootFinderBase, IRootFinder
    {
        #region Local variables and properties

        private Double d;
        private bool mflag = true;

        #endregion

        #region Interface variables and properties

        /// <summary>Gets a value that specifies whether the current function setup is solveable with the selected root-finding method.</summary>
        public override bool IsSolvable
        {
            get
            {
                return !(this.Inputs.fxnValInit * this.Inputs.fxnValPrev1 > 0);
            }
        }
        /// <summary>Gets and sets the inputs required for the root-finding algorithm.</summary>
        public override RootFinderInputs Inputs
        {
            get
            {
                return base.Inputs;
            }
            set
            {
                base.Inputs = value;
                if (Math.Abs(base.Inputs.fxnValInit) < Math.Abs(base.Inputs.fxnValPrev1))
                {
                    swap(ref base.Inputs.rootInit, ref base.Inputs.rootPrev1);
                    swap(ref base.Inputs.fxnValInit, ref base.Inputs.fxnValPrev1);
                }
                base.Inputs.rootPrev2 = base.Inputs.rootInit;
                base.Inputs.fxnValPrev2 = base.Inputs.fxnValInit;
                this.mflag = true;
            }
        }

        #endregion
        #region Interface methods

        // Local methods
        private void swap(ref Double val1, ref Double val2)
        {
            Double newval = val2;
            val2 = val1;
            val1 = newval;
        }
        /// <summary>Gets the root of the specified function.</summary>
        /// <returns>Returns the root of the function.</returns>
        public override Double GetRoot()
        {
            if (!this.IsSolvable)
                throw new Exception("This root-finding method requires an interval that contains the function root value.");
            if (this.Inputs.fxnValInit == 0.0) return this.Inputs.rootInit;
            if (this.Inputs.fxnValPrev1 == 0.0) return this.Inputs.rootPrev1;

            for (Int32 i = 0; i < this.MaxIterations; i++)
            {
                this.RootEstimate = this.EstimateRoot();
                this.FunctionValue = this.Function.GetValue(this.RootEstimate);
                if (this.IsConverged()) return this.RootEstimate;
                this.UpdateInfo();
            }
            throw new Exception("The root-finding method failed to achieve a solution.");
        }
        /// <summary>Estimates the root. This is used during iterations to calculate closer and closer estimates of the root.</summary>
        /// <returns>Returns an estimate of the root. This is used during iterations to calculate closer and closer estimates of the root.</returns>
        public override Double EstimateRoot()
        {
            Double s;
            if (this.Inputs.fxnValInit != this.Inputs.fxnValPrev2 && this.Inputs.fxnValPrev1 != this.Inputs.fxnValPrev2)
                s = this.Inputs.rootInit * this.Inputs.fxnValPrev1 * this.Inputs.fxnValPrev2 / ((this.Inputs.fxnValInit - this.Inputs.fxnValPrev1) * (this.Inputs.fxnValInit - this.Inputs.fxnValPrev2))
                    + this.Inputs.rootPrev1 * this.Inputs.fxnValInit * this.Inputs.fxnValPrev2 / ((this.Inputs.fxnValPrev1 - this.Inputs.fxnValInit) * (this.Inputs.fxnValPrev1 - this.Inputs.fxnValPrev2))
                    + this.Inputs.rootPrev2 * this.Inputs.fxnValPrev1 * this.Inputs.fxnValInit / ((this.Inputs.fxnValPrev2 - this.Inputs.fxnValPrev1) * (this.Inputs.fxnValPrev2 - this.Inputs.fxnValInit));
            else
                s = this.Inputs.rootPrev1 - this.Inputs.fxnValPrev1 * (this.Inputs.rootPrev1 - this.Inputs.rootInit) / (this.Inputs.fxnValPrev1 - this.Inputs.fxnValInit);

            // Use bisection if 1 of 5 conditions are met...
            Double beg = (3.0 * this.Inputs.rootInit + this.Inputs.rootPrev1) / 4.0;
            Double end = this.Inputs.rootPrev1;
            if (beg > end) this.swap(ref beg, ref end);
            if (!(beg <= s && s <= end)
                || (this.mflag && (Math.Abs(s - this.Inputs.rootPrev1) >= Math.Abs(this.Inputs.rootPrev1 - this.Inputs.rootPrev2) / 2.0 || Math.Abs(this.Inputs.rootPrev1 - this.Inputs.rootPrev2) < this.Tolerance))
                || (!this.mflag && (Math.Abs(s - this.Inputs.rootPrev1) >= Math.Abs(this.Inputs.rootPrev2 - this.d) / 2.0 || Math.Abs(this.Inputs.rootPrev2 - this.d) < this.Tolerance)))
            {
                s = (this.Inputs.rootInit + this.Inputs.rootPrev1) / 2.0;
                this.mflag = true;
            }
            else
                this.mflag = false;

            return s;
        }
        /// <summary>Determines whether a root for the function has been found.</summary>
        /// <returns>Returns true if the root for the function has been found; otherwise, returns false.</returns>
        public override bool IsConverged()
        {
            return (Math.Abs(this.FunctionValue) < this.Tolerance || Math.Abs(this.Inputs.rootPrev1 - this.Inputs.rootInit) < this.Tolerance);
        }
        /// <summary>Updates additional information used for solving the root-finding problem.</summary>
        public override void UpdateInfo()
        {
            this.d = this.Inputs.rootPrev2;
            this.Inputs.rootPrev2 = this.Inputs.rootPrev1;
            this.Inputs.fxnValPrev2 = this.Inputs.fxnValPrev1;

            // Assign to the proper side
            if (this.FunctionValue * this.Inputs.fxnValInit >= 0)
            {
                this.Inputs.rootInit = this.RootEstimate;
                this.Inputs.fxnValInit = this.FunctionValue;
            }
            else
            {
                this.Inputs.rootPrev1 = this.RootEstimate;
                this.Inputs.fxnValPrev1 = this.FunctionValue;
            }

            // Ensure b is the best root estimate
            if (Math.Abs(this.Inputs.fxnValInit) < Math.Abs(this.Inputs.fxnValPrev1))
            {
                swap(ref this.Inputs.rootInit, ref this.Inputs.rootPrev1);
                swap(ref this.Inputs.fxnValInit, ref this.Inputs.fxnValPrev1);
            }
            this.RootEstimate = this.Inputs.rootPrev1;
            this.FunctionValue = this.Inputs.fxnValPrev1;
        }

        #endregion
    }

    #endregion
}


